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Abstract We begin the construction and the analysis of nonosdllatory shock capturing 
methods far the approximation of hyperbolic conservation laws. These schemes 
share many desirable prop er ties with total variation diminishing schemes, 
but TVD schemes have at most first order accuracy, in the sense of 
truncation error, at e x t rem a of the solution. In this paper we construct 
a uniformly second order approximation, which is nonosdllatory in the sense 
that the number of e xtrem a of the discrete solution is not increasing in time. 

This is achieved via a nonosdllatory piecewise linear reconstruction of the 
solution from its cell averages, time evolution through an approximate 
solution of the resulting initial value problem, and averaging of this 
approximate solution over each cell. 


AMS-MOS Classification: Primary 65 MIO, Secondary 65 MQ5 
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1. Introduction. In this paper we consider numerical approximations to weak 
solutions of the scalar initial value problem (IVP) 

(1.1a) u t + f(u) x = u, + a(u) u x = 0 

(1.1b) u(x,0) = uo(x) . 

The initial data iio(x) are assumed to be piecewise-smooth functions that 
are either periodic or of compact support. 

Let 3 v h (xj, t„), xj » jh, t n = nr, denote a numerical approximation in 
conservation form 

(1.2a) vf*' - v} - - fj-yi) m (E„ ■ v"). 

Here E h is the numerical solution operator; X = r/h, and fj+yz, the 
numerical flux is a function of 2k variables 

(l-2b) fj+v 2 = fiyj-k* i,—.v"+*) 

which is consistent with (1.1a) in the sense that 
(1.2a) /(«,u,...,u) = /(u) . 

We consider the numerical approximation v h (x,t) in (1.2) to be a 
piecewise-constant function 

(1.3) v A (x,r) “ v/. x )-\n<*< X j+V2* nr < t s (n + 1)t . 

Accordingly we define its total variation in x to be 

(1.4) TV(v”) = = 2 ! v 7+i - . 
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If the total variation of the numerical solution is uniformly bounded in h for 

0 srsr 

(1.5) mv*(v)) * C • iT(uo) 

then any refinement sequence h - 0, r = o(A) has a subsequence fy - 0 so 
that 

(1.6) v hj —> u 
where u is a weak solution of (1.1). 

If all limit solutions (1.6) of the numerical solution (1.2) satisfy an entropy 
condition that implies uniqueness of the IVP (1.1). then the numerical scheme is 
convergent (see e.g [3] y [12]). 

Recently we have introduced the notion of total variation diminishing (TVD) 
schemes (see [3]), where the approximate solution operator is required to dimin¬ 
ish the total variation (1.4) of the numerical solution at each time-step 

(1.7) 7V(v J,+1 ) <s TV(v H ) ; 

these schemes trivially satisfy (1.5) with C = 1. Some early work along these 
lines was done by van Leer in [15]. 

TVD schemes are non-osdllatory in the sense that the number of local 
extrema in the numerical solution is diminishing in time (as is customary we use 
"diminishing" loosely as short for "non-increasing", throughout this paper). 
Moreover, the value of an isolated local maximum may only decrease in time, 
while that of a local minimum may only increase. 

We were able to construct TVD schemes that in the sense of local truncation 


error are high-order accurate everywhere except at local extrema where they 
necessarily degenerate into first-order accuracy (see [4], [13], [10], [11], [14]). 
The perpetual damping of local extrem a determines the cumulative global error of 
the "high-order TVD schemes" to be 0(h) in the L» norm, 0(/i ;v2 ) in the L 2 
norm and 0(h 2 ) in the L\ norm (see [17]). 

In this paper we introduce a larger class of non-osdllatory schemes, in which 
the solution operator is only required to diminish the number of local e xtrem a in 
the n umerical solution. Unlike TVD schemes, which are a subset of this class, 
non-osdllatory schemes are not required to damp the values of each local 
extremum at every single time-step, but are allowed to occasionally accentuate a 
local e x tr em u m . 

hi a sequence of papers, of which the present paper is the ilrst, we show how 
to construct non-osdllatory schemes that are uniformly high-order accurate (in 
the sense of global error for smooth solutions of (1.1)). In this first paper we 
describe a second-order accurate scheme of this type. 

The fact that the number of local ex tr em a in the numerical solution may only 
diminish in time is sufficient by itself to guarantee that the application of the 
scheme to monotone data results in a monotone function. Thus non-osdllatory 
schemes, like TVD schemes, are monotonidty preserving. In particular, when 
applied to a step-function, they do not generate spurious oscillations. 

We note that since the number of local extrema in the solution of non- 
osdllatory schemes is bounded by that of the initial data, uniform boundedness of 
its total variation (1.5) follows immediately if the maximum norm of the solution 
is shown to be uniformly bounded. 
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2. Design Principle and Overflew 

In this section we describe how to construct a non-osdllatory scheme that is 
uniformly second-order accurate. 

Integrating the partial differential equation (1.1a) over the computational cell 

(*/—V2» */+i/ 2 ) x ('»> Vt-i) we get 

( 2 . 1 a) « 7 + ‘ “ H) ~ ») " fj-l/M] 

where 

( 2 . 1 b) }j+vM - 7 0 ) * » 

and 

( 2 - lc ) “ 7=7 “(*.'*) * • 

We observe that although (2.1a) is a relation between the cell-averages uj 
and uf+ l , the evaluation of the fluxes in ( 2 . 1 b) requires knowledge 

of the solution itself and not its cell-averages. 

As in Godunov’s scheme and its second-order extension by van Leer [16] and 
Colella and Woodward [2], we derive our scheme as a direct approximation to 

(2.1) . We denote by vj the numerical approximation to the cell-averages uf of 
the exact solution in ( 2 . 1 c) y and set vf to be the cell-averages of the initial data. 
Given v" = {vf) we compute v " +1 as follows: 

First we reconstruct u(x,/„) out of its approximate cell-averages {vf) to the 
appropriate accuracy and denote the result by I(x; v"). Next we solve the IVP 

( 2 . 2 ) v, + /(v) x = 0 , v(x, 0 ) = L(x; v") 
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and denote its solution by v(x,r). Finally we obtain vj l+1 bv taking cell- 
averages of v(x,t) 

(2.3) v?+- = | v(i,t) dr . 

The averaging operator in (2.3) is non-oadllatory, therefore the number of 
local extrema in v rt+1 (interpreted as a mesh-function or the piecewise-constant 
function (1.3)) does not exceed that of v(x,t). Assuming v(x,t) to oe the exact 
solution of (2.2) implies (since the exact solution operator is TVD) that the 
number of local extrema in v(x,t) is less than or equal to that of 
v(x,0) = L(x; v*). Therefore if the number of local extrema in L(x; v n ) does 
not exceed that of v*, then the resulting scheme is non-oscillatory. 

We conclude that the design of non-oscillatory high order accurate schemes 
essentially boils down to a problem on the,level of approximation of functions: 
Given cell-averages uj of a piecewise-smooth function u(x), reconstruct u(x) 
to a desired accuracy. Prior to studying this problem we tackle another related 
question in approximation of functions, that of constructing a non-oscillatory 
high-order accurate interpolation of piecewise-smooth functions. 

In section 3 we construct a non-oscillatory piecewise-parabolic function 
Q(x; u ) that interpolates a piecewise-smooth function u(x) at the mesh points 

(2.4a) Q(xj; u) = u(xj) 

and satisfies, wherever u(x) is smooth, 

(2.4b) Q(x\ u) = u(x) + 0(h 3 ) 

(2.4c) ■— C(x ± 0; u) = — m(x) + 0{h*) . 

- 7 - 




V 






In section 4 we make use of this non-osdllatory piecewise-parabolic interpo- 
lant io design a non-osdllatory reconstruction of a piecewise - s mooth function 
from its cell-averages. As in [16], [2], [5], ind [9] we take L(x\ u) to be the fol¬ 
lowing piecewise-linear function 

(2.5a) L(x; u) = uj + Sj(x - xj)Jh for \x - xj\< h/2 , 

Unlike the above references that present "second-order accurate" TVD 
schemes, we compute the slopes SJh from Q(x; u) by 

(2.5b) S/h = Q(xj - 0; u), Q(xj + 0; u)j . 


Here m(x,y ) is the min mod funoion 
(2.6) m(x 


'■*) = jo 


s • min(|x|,H) if sgn(x) = sitn(y) 
. otherwise 


We show in section 4 that L(x; u) is a proper reconstruction of u(x) in the 
sense that whenever u(x) is smooth 

(2.7a) L(x; 5) = «(*) + Off. 2 ) 

and 


(2.7b) 


L(x; u) = u(x) + 0(A 3 ) . 


Here u(x) = h~ l u(x ■+■ y) dy and L(x; u) = 

* h~ l L(x + y; u) dy; like Q(x; u), the latter is also a non-osdllatory 
piecewise-parabolic interpolant of u(x), 

(2.7c) L(xj; 5) = u(xj) . 
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We remark that the "second-order accurate" TVD schemes described in the 
above mentioned reference use a slope Sj/h in (2.5a) that approximates 
(d/dx) u(xj) to 0(h ), and their loss of second-order accuracy at local e xt rem a 
points is due to lack of smoothness of the coefficient in the 0(h) term at these 
points 1 . This problem is circumvented in the present scheme by taking Sj/h to 
be (2.5b) which is an 0(h*) approximation to (d/dx) u(xj). Unfortunately there 
is a price to pay for this extra accuracy, namely the loss of the TVD property. As 
in TVD schemes 

(2.8) JV(v« + >)s7V(i(;v')), 

however here 

TV(L('\ v")) a TV(v n ) 

and indeed the scheme may occasionally increase the variation of the numerical 
solution. Although we prove that the scheme is non-osdllatory we have not been 
able as yet to complete a proof of uniform boundedness of the total variation of 
the numerical solution; this is due to lack of techniques to verify uniform bound¬ 
edness of the maximum norm of the numerical solution. 

In section 5 we study the proposed scheme in the constant coefficient case. 

We verify that it is uniformly second-order accurate, examine its behavior at local 
extrema points and get estimates for the possible increase in total variation per 
time-step. 

In this paper where we consider numerical schemes of the form (1.2) that ar- 


1. We repeat that the results of [8] and [11] imply that TVD schemes, no matter how 
they are constructed, must have this loss of accuracy at local extrema 


derived from approximating the relation (2.1), it is only natural to consider trun¬ 
cation error in the sense of cell-averages, i.e. we say that the scheme (1.2) is 
second-order accurate if 

(2.9) S"* 1 - E» -5" + 0(A 3 ) 

where 5" is the cell-average (2.1c) of the exact solution. Since 

(2.10) u(x) - u(x) + 0(h~) 

whenever u(x) is smooth, (2.9) holds also for pointwise values of the solution. 
However, in the cont e x t of 3rd and higher order accurate schemes, this difference 
in definitions of truncation error will be not only conceptual but of practical 
importance as well. 

Up to this point we have assumed that v(x,t) in (2.3) is the exact solution 
to (2.2). The resulting scheme 

(2.11a) vj* 1 = v}~ klfj+^v) -/,-„&)} , 

where fj+ i/:( v ) i» (2.1b) applied to »(*,/), 

( 2 - llb ) /y+i/2<v) - “ /(v(x,/)) dt , 

is certainly second-order accurate in the sense of (2.9). Starting with the exact 
cell-averages vj - uj in (2.11) we get from (2.7a) that 

(2.12a) v(x,/) ” a(x,f + f*) + 0(h*) for OsfSr 

and consequently 

(2.12b) - }j+vM + 0(h *>. 

which implies (2.9) due to the sufficient smoothness of the coefficient in the 
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Oih 2 ) term in (2.12b). 

1 a section 6 we replace the exact solution v(x,f) in (2.3) by an approximate 
one, which we denote by v„(x,/) . This approximate solution is conservative. 
TVD, and second-order accurate in the sense of (2.12a). Thus replacing v(x,f) 
in (2.3) by this approximate solution results in a conservative scheme that is non- 
osdllatory and uniformly second order accurate. 

We remark that an alternative approach to the above is to approximate 
/;+ inH v ) ^ (2.11b) by using a midpoint rule (or trapezoidal rule) for the 
integral and by replacing v(x,f) with a non-oscillating second-order accurate 
approximate one v„(x,r) (see [16] and [2]). The resulting scheme 


(2.13a) 

vj*' = 

(2.13b) 

}j*Vl ” /Wf/.l/I. T/2)) 


is certainly second-order accurate, and it is -non-oscillatory in the constant coeffi¬ 
cient case. Since we have not used the cell-averaging (2.3) to derive this scheme, 
we cannot ascertain in general that the resulting scheme is non-osdllatory. 
Nevertheless, our numerical experiments as well as many other experiments in 
the context of TVD schemes (see e.g. [1], [2]) demonstrate that the numerical 
results are non-osdllatory in many (if not all) applications. 

In section 7 we present some numerical e xp e r i ments that compare the present 
scheme with a typical "second-order accurate" TVD scheme. 

3. Nonoadllatory Interpolation. 

The oscillatory nature of second order accurate Lax-Wendroff type schemes 
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results from a Gibbs phenomenon associated with high-order interpolation across 
discontinuities. In this section, as a preparatory step towards Signing a nonos- 
dllatory approximation to (1.1.), we construct a non-oscillatory piecewise- 
parabolic interpoiant Q(x; u) to a piecewise-smooth function u(x) such that 

(3.1a) Q(x,; u) = u(x,) 

(3.1b) Q(x,u) m g i+u2 (x;u), x { * x s x (+1 , 

where ft+ 1/2 i* a quadratic polynomial, and 

(3.1c) fi(x; u) - u(x) - O^ 3 ), -j- Q(x ± 0; 11) - u(x) - Oih 2 ) 
whoever u(x) is smooth. 

Q(x; u) is non-osdllatory in the sense that the number of its local e xtr em a 
does not e x ceed that of u(x). 

Since 

qi+i/2(*<; “) = “i* li+vito+ii u ) y tt i+i 
it can be written in the form 

(3.2a)q / +i /2 (x;u) « u, + d ( +i/2 u-(x - x t )/h + J0/+1/2 u (x - x,)(x - Xt+J/h 2 

where 

(3.2b) di+vi “ = “i+i “ u i 

and D iJr]f 2 u is yet to be determined. 

0.+ 1/2 u = * 2 li+yifa “). x, ^x*x, + i . 
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We consider as candidates for < 7 ^ 1/2 the two quadratic polynomials g, 
and qi+ lt interpolating u(x) at fo-j, x h x i+1 ) and (x h x ;+1> x l+ 2 ), respec¬ 
tively, and choose q t + ^ to be the one that is least oscillatory in [x,-, x (+1 ]. 
Both qj , j = i and j ~ i + 1, can be written as (3.2a) with u = D j u 

where 


(3.2c) Dju = dy + 1 / 2 j< - - 2uj + . 

Since the least oscillatory of q t and can be characterized as the one that 
deviates the least from the line connecting (x^u,) with (xj+i,Uf+i) we choose 
u “ (3.2a) to be 

(3.2d) D i+U2 u = m(D,u, D i+ 1 u) , 

where m(x,y) is the min mod function 


(3.3) 


m(x,y) 


s • min([x|,lyD .if sgn(x) = sgn(y) = s 
0 ' otherwise. 


If u(x) is smooth in [xy_ lt xj* 1 ], then qj as a quadratic interpolant of u 
satisfies 


(3-4)^(x) - u(x) = 0(A 3 ), -£■ qj(x) - ■£■ u(x) = Oih 2 ), Xj - x s: x s: x >+1 . 

If D,u • D j.\u i 0 then $ 1 + 1/2 is either $j or Otherwise we set 

A+ 1/2 x = 0, but then smoothness of u implies that Dju = 0(h 3 ) and conse¬ 
quently qi+y2 “ qj ~ 0(h 3 ) for j = i, i + 1. Thus (3.1c) follows from (3.4). 

We turn now to prove that Q(x\ u) is a nonosdllatory interpolant of u, 
i.e. that the number of its local extrema does not exceed that of u. We do so by 
showing a oqe-to-one correspondence between local e x t rem a of Q to those of the 
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mesh function {uy}, the number of which certainly does not exceed that of u(x). 

Q may have a local extremum in either the interior of some interval 
(x,,x l+1 ) or at a mesh point x,. The first case, which will be refered to as 
interior-extremum, occurs when there is a point x*, x, < x* < Xj+ lf such that 

Q(x*‘, «0 = 0, but -£■ < 7 1>1/2 f 0 . 

From (3.2a) if follows that Q has an u tenor-extremum. m (x„ x^i) if and 
only if 

(3.5) |0J+l/2 “I > 2 l4-»-l/2 m ! • 


qi+i /2 = ?,+i/ 2 CO» toe valuc °f toe interior-extremum is then 


(3.6) 


1 - f 4+1/2“ lV 
* +w ~ «f - J Di+v* ' 2 


if Di+y 2 u < 0 it is a local maximum; if Di+yi “ > 0 it is a local minimum. 
Since Di+yju = m(D,u, D {+1 u), (3.5) holds if and only if 

(3.7a) Di u ' D f +in > 0 


(3.7b) \Pj u| > 2| di+y^ uj, j = i, i + 1 . 

This implies that qi+yi has a locai extremum in (x f , x i+1 ) if and only if both 9) 
and 9 ) 4.1 also have a local extremum in (x,, x 1+1 ) and of the same kind. Since 
a parabola has at most one local extrem um, it follows then that 9) does not have 
a local extremum in (Xf-j, x,) and 4+1 does not have one in (x,+i, x, +2 ). 
Consequently Q is monotone in both (xj« lt x,) and (x <+1 , x^^, but in an 
opposite sense, i.e. 4 - 1/2 M ‘ 4 + 3/2 M < 0 ; toe latter implies that u has a local 
e xtr e m um in [x,-, x,+i] and that either u,- or u ( 4.1 is a local extremum of the 
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mesh function {uj} (for obvious reasons the case u, = u^i is counted as a 
sipglg -^TiT ffmiim ). The above analysis also shows that In tenor-extrema are iso¬ 
lated, i.e. if Q has an interior-extremum in (x <t x (+1 ), then it is the oniy local 
extremum of Q in (x,-!, x l+2 )- 

We turn now to examine the case that Q has a local extre m u m at a mesh 
point x,; this will be refered to as a mesh-extremum. The above observation 
t hat interior extrema are isolated excludes the possibility that Q has an interior- 
extremum in either (x,-j, x,) or (x,, x, + i) and consequently Q is monotone in 
tire* intervals. This im plies that u * di+ 1/2 a < 0 and therefore is a 

local extr emum of the mesh function {u^}. This concludes the proof that Q(x\ u) 
is a non-oscillatory interpolant of u. 

We next express the non-osdllatory nature of Q in terms of total variation. 
If |D;+i /2 M l ^ 2|dj+ 1/2 “! *en (2.5) implies that Q is monotone in (x ; , xj+J. 
Thus 


(3.8a) \Dj+in “I 25 2 M/-m /2 “I =* TV\x,^(Q) ~ ty+1/2 “I - 

If |D, + i/2 M i > 2 M/+i/2 “I then Q has a local extremum in (x Jt xj+i) and 

TVfo^CQ) = \q]+v 2 ~ u j\ + i“;*i ~ 4j^\n\ • 

Using (2.6) we get 

(3.8b) 1 Dj+y2 u\ > 2\d<+ m “I => TV [Xj Jj ^(Q) 


= |d/+i/2 “I + i^+l/2 M l 


u 

d J+1J2 u 



We conclude that 


(3.8c) 0 s 7V(fi) - 2 1/2 “I * 2 lD—l/2 “I f - i) 

j miM + 2J 

S T 2 l^m+1/2 “i • 

The sum in the RHS of (3.8c) is taken over the set of indices M of intervals 
(xm, **+i) in which |D m+1/2 “I > 2K.+V2 «l» where Q has interior- 
extrema. 

Next we show that if u(x) is a piecewise-smooth function of bounded varia¬ 
tion, then 

(3.9) lim 7V(Q( ;u)) = 7V(u) . 

A-0 

We observe that in this case the number of intervals in M is finite and is uni¬ 
formly bounded by the number of local extrema in u(x). Hence (3.9) will follow 
if we prove that 1/2 tt 0 as A - 0 for all m £ M. To accomplish this we 
show that for h sufficiently small M does not include intervals {x m , in 

which u(x) is discontinuous. To see that let us examine the case where u(x) 
has a discontinuity at x € (x h x 1+1 ) . Qearly d i+1/2 u approaches the size of 
the jump in u while y^u approaches zero as /i - 0. Consequently 

(3.10a) p i u/d i ^ y2 “I = |1 - d i-\n* ld i+\n “1-1 ** * - 0 
Hence for h sufficiently small 
(3.10b) 2|dj +1/ 2 m| > pi u| Pi+vi u| 

which implies i g M. 
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4 . Non-oadllatory Reconstruction. 

Let u(x) be a piecewise-smooth function and denote by u(x ) its mean over 
(x — h/2jc + A/2), i.e. 

(4.1) “to « j “to dy . 

We denote uj = u(xj) and refer to these values as cell-averages of u(x). Given 
{£)}, the task in hand is to reconstruct u(x) to 0(h?) in a non-oscillatory way; 
denote the approximately reconstructed function by L(x; u). To achieve 0(h?) 
accuracy it is sufficient to consider L(x; u) to be a piecewise-linear function. To 
make L(x\ u) a non-osdllatory approximation we use the non-osdllatory piece- 
wise parabolic interpolation Q(x; u) to compute its slopes as follows: 

(4.2a) L(x; u) = uj + S)(x - xj)Jh for \x - xj\ < y A 

(4.2b) Sj = h- Q(Xj - 0; i), Q(x t + 0; S)| . 

Here m is the min mod function (3.3); d i+ y 2 u and u are (3.2b) and 

(3.2d), respectively. 

We note that L(x; u) may be discontinuous at {xy +1 ^} and that 
(4.3a) L(xj; u) = uj . 

To see that wherever u(x) is smooth 
(4.3b) L(x,u) - u(x) = Oik 2 ) 

we observe that in this case 

(4.4a) u(x) = u(x) + (^(A 2 ) 
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and therefore it follows from (3.1c) that 


(4.4b) J S l = ~Zi + 0(A2) “ * + 0<h2) 

Consequently the RHS of (4.2a) can be expanded as 
(4.4c) L(x; u) = uj + (x - xj) u(xj) + O^h 2 ) 

= u(x) + (X* 2 ) for lx - xj\ < y h 
and thus (4.3b) follows. 

Denote by L(x; u) the mean value of L(x; u) in (x - hi2, x + h/2), i.e. 
(4.5a) L(x\ S) - j L(y\S)dy. 

Using (4.2a) to evaluate the integral in (4.5a) we find 

(4.5b) L(x; u) = uj + dj+w u 

• (x - xj)/h + (V2)(Sj+i - Sj)(x - xj)(x - xj+^h 2 , 
for xj ^ x s xj+i 

(4.5c) L(xj\ u) = uj . 

Hence L(x; u), like Q(x; u), is a piecewise-parabolic interpolant of u(x). 
Comparing (4.5b) with (3.2) we find that for x^sxs xj+i 

(4.6a) L(x; u) - G(x; u) = y(5/.-ri “ Sj - Dj+ia **)(* “ x j)( x “ */+iy* 2 • 
From (4.4b) we see that 5, = h u(x,) + 0(h 3 ) (Note that this is true 
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also at local extrema points) and therefore 

S J+ 1 " S J = hZ “(*/+i/2) + 0(h 3 ) . 
On the other hand (3.2) shows that 

-> d 2 

D )+\n u = h 2 **(*./+1/2) + 0 (A 3 ). 

Therefore 


^ " Sj - D y+y2 u = 0(A 3 ) 

which shows that RHSof (4.<5a) is 0(h\ Since (2.1c) show, that 

Q(x; u) - u(x) = 0(h 3 ) 
we conclude from (4.6a)-(4.6b) that 

(4,6c) fa; •*) “ u(x) = 0(h 3 ) . 

We turn now to prove that £(*; S) is a non-osdllatory approximation to 
J(x); this certainly implies that £(x; 5) is a non-o*ffla,ory approximation to 
»(*). We shall do so by showing that TV^jdfr; i)). the totai-variauon of 
L ( x * **) [*/> x ]+ iJ, which has the value 


(4.7a) i)) - ^ (£y| + + |rfy +1/2 u - 

can also be expressed as 


y(Sy + Sy + i)l • 


(4.7b) 


7 V k^-j( £ ('; “» - “«(K+w <fl, so. 


Then it follows immediately form (4.7b), (4.3a) and (3.8) that £ is monotone in 
[xy.*, + 1 ] if and only if Q is; consequently £ is a non-osdllatory approximation 
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to u(x) in exactly the same sense as Q is to the interpolated function (see sec¬ 
tion 3). 

Next let us denote 

(4.8a) Sf = h ■ Q(xj ± 0; i) , 

i.e. 

(4.8b) Sf = dj-y 2 u + y Dj-yi u. Sf = d J+ ^ u ~ y Dj+yz u , 
and observe that (4.2b) implies that 

(4.8c) j(|S,| + |S y+1 |) = Sf) | + |m(5^ lf S/+i)|] 

^ yOsf I + |S/+iD = y |m./+i /2 “ - y Dj+m. “1 + \ d j*vi u + y Dj+vi « 
= max||rfj4.i/2 S|, ylD ;+1/2 ulj • 

We note that if ld/ +1 /2 **| ^ l^l^j+i/2 “1 then 

sgn|^i /2 « ± y D j+in “j * sgn(^/+u 2 “) a 0 

which in turn implies 

sgn(5j) * sgn(d)+i/2 u) 2* 0, sgn(.Sy +1 ) • sgn(d) +1/2 u) ^ 0 
It follows then from (4.8c) that the RHS of (4.7a) is \dj+y2 5]. This shows that 

(4.9a) | dj+y2 “1 ^ y )Pj+V2 “1 =* TV [“)) = M/+1/2 “I * 
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To complete the verification of (4.7b) we still have to show that 


(4.9b) I d J+M S| < i 1|D,+ W si = TVj^jCK-; 5)) = y |D, +1 5] . 

First we observe that 

(4.10) S+ - Si~ = (di+ifi “ “ y ^i+1/2 “) “ (di—i /2 “ + Y D i-i/l “) • 

= A “ ” Y (D ^^ “ + **) 

Since (3.2d) implies 

(“•ID |O(S|*Y(|Oi-w5| + |Di + w50 

we conclude from (4.10) that 

(4.12) (5, + - SD ' HJpPi «*) 0 • 

We turn now to prove (4.9b). First let us consider the case that Q(x; u) 
has a local maximum in (xj, xy +1 ), i.e. Dj u < 0, Dj +1 u < 0, and 

\4j+ 1/2 “1 < Y lty+1/2 “I* 

It follows from (4.12) that 

(4.13a) Sj~ a Sf = d)+y 2 U - y Dj+yi u > 0 

(4.13b) 0 > u + y ty+i/2 “ = ty+i ^ Sj+i • 

The relations (4.13) and the definitions (4.8a), (4.2b) imply that 
(4.14a) Sj = S* = dj+M “ “ y 0y+i/2 u 
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V 


(4.14b) 


S J +1 * S J~+1 = d j+\n “ + 2 D J+vi “ • 

The same analysis shows that (4.14) holds also for the case that Q(x; u) has a 
local minimum in (x ; , (4.9b) follows immediately from (4.14) and (4.7a). 

We note that since L(x; u) is continuous at xj 

(4.15) 7V(L( ;u» = 2 u)) = 2 max^ +w iTl, | \D J+l/2 S\ 

- J My+i/i “1 + Pm+vi "I " Ki+yj “ij • 

Here M is the set of indices of intervals (x^^.^) in the interior of which L 
(and also Q) has a local extrem um. The number of these intervals is finite and 
is bounded by the number of local e xtrem a of u(x). Comparing (4.9) with (3.8) 
we note that 

(4.16) 7V(L(;5))a7y(e(;5)). 

5. The constant coefficient case. 

In this section we study the constant coefficient case 
(5.1) Uf + <zUj = 0, a = const. 

The exact solution of the IVP (2.2) is 

(5-2) v(x,r) = L(x - at ; v ") . 

Hence our scheme (2.3) is 

(5.3a) v/ +1 = -jj- JL(x - at; v") dx = L(xj - ar; v") 
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where L is (4.5a). We have shown in section 3 that the number of local 
extrema in L(x\ v n ) does not exceed that of v n . Since v fl ' rl in (5.3a) is a cell* 
average of L, it follows that the number of local ertrema in v ** 1 does not 
exceed that of v", and consequently the scheme (5.3a) is non-oscillatory. 

Using (4.5b) in (5.3a) we get the following expression for the scheme 

(5.3b) v/ +1 * (£» • v")y 

yf - pdj-yi v” - 1/2 |i(l - (0(57 - 57- 1 ) 
vj - pdj+vi v" + 1/2 |i(l + |i)(5”+i - Sf) 

E h denotes the operator form of the finite difference scheme; 

CFL-number, is assumed to satisfy 

(5.3c) W*l. 

We turn now to prove that (5.3) is second order accurate in the sense of 
(2.9), i.e. if u"(x) denotes the mean value (4.1) of u(x,/„) then 

(5.4a) 57* 1 - (£* • i!")/ - 0(A 3 ) . 

To show that we observe that in the constant coefficient case (5.1) 

= iT(xj - at), and by (5.3a) (E h • fV T ) j = L(xj - aT; u"). Hence the LHS 
of (5.4a) is nothing but 

(5.4b) u'ixj - at) - L(xj - ar; iT) , 

which is 0(h 3 ) as a direct consequence of (4.bc). 

Next we study the time-dependence of the total variation and the maximum 
norm of the numerical solution (5.3). In section 2 we have pointed out that 


if a > 0 
if a<0 ; 

U- = Xa, the 
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(5.5a) 7V(v" +1 ) * TV(L(*; v")) . 

Using (4.15) and (5.5a) we get the following upper bound on the possible growth 
of the total variation of the numerical solution per time-step 

(5.5b) 7V(v» + ‘) - 7V(v") * | D m + M v»| - K. +1/J v»|| . 

Here M„ is the set of indices of intervals (x m , in the interior of 

which L(x\ v n ) (and also Q(x; v")) has a local extremum. The number of these 
intervals is finite and remains uniformly bounded in time by the number of local 
extrema in the initial data. 

Clearly the upper bound (5.5b) is overly pessimistic. It estimates the possi¬ 
ble increase in variation in the reconstruction step due to replacing the cell- 
averages vj by the piecewise-linear function L(x\ v*). It does not talce into 
account the possible decrease in variation in the averaging step (2.3), resulting 
from doing just the opposite, i.e. replacing the piecewise-linear function 
L(x - or; v n ) in (5.2) by its cell-averages (5.3a). 

In the following we shall examine the temporal behaviour of the local 
extrema of the numerical solution and its total variation by analysing the explicit 
values of the cell-averages given by (5.3b). To simplify our presentation 
let us assume that a > 0. 

First we note that (4.8b) implies 

(5.6.) \Sj ~ Vll s W + IVll s 2 ““(14-l/J ""I. \ |Dj-u2v"|) . 

Hence 
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(5.6b) v"| a - 1|D,- W v”! =, hj-vi I - I®/ “ Sf-il/ty-i/i v"| s 2 

Rewriting (5.3) in this case as 

(5.7a) v ;- 1 - v; - 4 . 4_ yJ v” - V2 4.0 - nty-yi 4-yi v" 

- (1 - < 4 - 1 / 1 ) «; * <»;-!/2 v ”-l 

where 

(5.7b) <*y-i /2 ~ M- + y M-(l - M-by- 1/2 • 

Wc see that the CFL condition 0 < u ^ 1 and (5.6b) imply that 
(5.7c) 

thus we conclude 

(5.8) 14-1/1 v"! ^ Y lO/.y, V*| => y7 +1 ( (v;-„ vjn . 

Relation (5.8) shows that if v" is monotone for J L s j ^ /*, i.e. 
vy t * v 7 t *i * • • • * v y-f or Vy t * vy tM - 2 v y- , then V-+1 is mono¬ 

tone for Ji + 1 i j £ /*, and in the same sense. Relation (5.8) also shows 
that mesh-extrema of v*, i.e. those for which Q has ici local extremum at a 
mesh point, are being damped at the n-th time-step. Namely, 

(5.9a) \d )±v2 v"\ 2 y \D J±U1 v-| f vJ.jSv/a: =» max^;* 1 , vf*; 1 ) s vj 

(5.9b) !d ;±1 /2 v"| 2 y lD y±L r 2 w"|, vf +1 =» min^"' 1 , v^, 1 ) 2 vf . 

We turn now to consider interior local extrema of v n , i.e. those for which 
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Q has its local extremum in the interior of some (x ( , x (+1 ). We recall that such 
a .j e xtiem u m is characterized by |d|+i /2 v'*| < 1/2 (D^i n. v "i that 5" and 
5f +1 in this case are given by (4.14); therefore - J" - v". From 

(5.3a) and (4.6a) we see that in general 

(5.10a)v|'+ 1 1 - Q(*<+i - aT;v*) - y p.(l - |fc)(f>r+i/2 ~ $Pm + • 

Hence 

(5.10b) | d l+vl v”| < y |D, +w v"| * vftf = CC*i+i - «t; v”) . 

Relation (5.10b) confirms the second order accuracy of the scheme at local 
extrema. Although it does not necessitate accentuation of the extremal values, as 
vftf in (5.10b) may still be in [vf, vf +1 ], it dors allow vfo 1 to deviate from 
thi* interval by as much as 

(5.10c) y |D,+ W **|(k+VJ v"|/|D,*v 2 v”| - 1/2)' 

Thus (5.10b) is the essential difference between the present scheme and the 
"second order" TVD schemes. 

A s imilar analysis, which we do not present here, shows that if vf is a 

mesh-extremum then j — i, i -t 1, relates to Q(xj ~ ar: v n ) in the fol¬ 

lowing way: 

(5.11a) v/* 1 x Q(xj - ar; v"), J * i, i + 1, if vf is a maximum 

(5.11b) v/ +1 s C(xy - aT; v’), ; = /,/ + 1, if vf is a minimum . 

From (5.9)-(5.11) we deduce the following relation between the total varia- 
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tion of the numerical solution and that of its piecewise-parabolic interpolant Q: 
(5.12) 7V({C(i y - ar,»")}) s rv(v” + ') s TV(Q (-; v")). 

The LHS of (5.12) is the total variation of the mesh function 
{Q(xj - ar, v 71 )}. Relation (5.12) suggests to consider an equivalent definition 
TV of the total variation of the numerical approximation of the form 

7V(v") = rv(C(-; v")) 

with the hope that the scheme (5.3) is TVD with respect to this modified defini¬ 
tion. Unfortunately our numerical experiments have shown that there are 
instances, although rather rare, that TV(y ') is increasing with n; the same is 
true for 7V(v n ) = 7V(L(-; v*)). 

As we have mentioned in the introduction, because of the nonosdllatory 
nature of the scheme, uniform total variation boundedness of the numerical solu¬ 
tion is implied by uniform boundedness of its maximum norm. If we follow a 
particular local maximum of the initial data we see from (5.9)-(5.10a) that it 
actually decreases most of the time, and whenever it does increase (5.10c) and 
(3.10) suggest that it does so by a "small amount" that vanishes with h -> 0. 

Since the initial data is only piecewise-smooth we have not been able as yet to 
rigorize these arguments. 

We remark that our numerical experiments clearly indicate that in a normal 
computational situation the maximum norm of the numerical solution is indeed 
uniformly bounded. We feel that our inability to prove this fact stems only from 
lack of theoretical tools to analyse pointwise regularity of the numerical solution. 
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6. The nonlinear case. In this section we describe an approximate solution 
v„(x,t) of [5] for the IVP (2.2) 


(6.1) v f + /(v) T = 0, v(x.O) = L(x; v") . 

This approximate solution is consistent with the conservation form of the equation 

(6.1) in the sense that the cell-averaging (2.3) results in a scheme in conservation 
form i.e. 

(6.2) v/ +1 = j J[ v„(x,t) dx = vf - X(/y +ly2 - /y-!^) 

where the numerical flux fj+\n is consistent with /(u) in the sense of (1.2c). 
Furth er mo re , the approximate solution operator is TVD 

(6.3) 7V(v„(-; 0) = rV(v n (-J 0)) = TV(L(-; v")) for 0 * * s t 

and thus by the reasoning presented in section 2, the resulting scheme (6.2) is 
non-osdllatory. 

We turn now to outline the derivation of this approximate solution. To sim¬ 
plify our presentation we ignore entropy considerations and refer the reader to 
future papers for details of appropriate modifications. We assign to the point 
Xj+y2 a characteristic speed that corre s ponds to the Rankine-Hugoniot speed 
aj+y 2 of the two neighbouring cell-averages vf and v? +1 




(6.4a) 


a t + U2 


v 7 * v '?*' 

V J+ 1 V J 

a(?f) if vj* = v/ +1 


and denote by a(x) the piecewise-linear interpolant of {aj +1/2 }, i.e. 


(6.4b) a(x) = ajf-i/2 + ( 5 }+ 1/2 - ay-1/2) * (x “ Xj.^h 


for Xy_ !/2 5X5 X^i/2 . 

The approximate solution v„(x f t) is defined by specifying its constancy 
along the approximate characteristics 

(6.5a) x(f) = x 0 + a(xo) • t 

i.e. 

(6.5b) v II (x(f), t) m v,,(xo,0) = L(x 0 ; v H ) . 

Using (6.5a) and (6.4b) to ex p res s x 0 in terms of x and t 
6.5c) xq(x,0 = X]—\n + h [x - x y - w (0)/[x y+1 ^(0 - x ; _v2(0] 


for x^— 1 / 2(0 < x < x^ +L7 (0 
we get from (6.5b) that 


(6.6a) v n (x,0 


* 

= L xy_ 


1/2 


+ h 


* ~ Xy—1/2(0 

^+1^(0 “ X/—1/2(0 ’ 


for Xy_^(0 < X < Xy* i/2(0 


where 

(6.6b) x,+ 1 / 2(0 = x i+ i/2 + r • a,.,, i/2 • 

Taking cell-averages of the approximate solution (6.6) we get the conservation 
form (6.2) 

( 6 . 7 a) v} +: = vj - k(/y+i/2 ~ /y-1/2) 
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with the numerical flux 


(6.7b) 


h 


+ 1/2 “ 


f(vf) + 1/2 aj+i/z (1 - X aj-y 2 ) ‘ 

f(v"+ 1) “ 1/2 a} i-1/2(1 + V 5 }-^ 3/2) ' £/-: 


if 5)^ i/2 ^ 0 
if <*}-rl/2 s 0 


where 


(6.7c) i, = 57/(1 + X(5) +1 /2 - S,- u 2 )] • 


Note that (6.7) is identical to (5.30) in the constant coefficient case. 

We turn now to prove that the scheme (6.7) is uniformly second-order accu¬ 
rate in the sense of (2.9). Starting with the exact cell-averages vf = uj in (6.7) 
this amounts to showing that 


(6.8a) /jf+ 1/2 ~ fjrvM + Oih 2 ) 

with a sufficiently smooth coefficient in the Oik 2 ) term; here fj+yi is the 
n umerical flux (6.7b) computed with the exact cell-averages, and fj+ 1 / 2 OO is 
(2.1b). We shall do so in two steps: first we shall show that 

(6.8b) fj+v2(u) = ^\f(L(xj+y 2 ; u")) +/(I(x 0 (x /+U 2, T);iT))] + Oih 2 ) , 


where x<j(xj +1 / 2 ,t) is (6.5c), and then we shall verify that 

(6.8c) ±[f(L(xj+y 2 \ O) + /(I(xo(x y+ i/ 2 , t); 5"))] = f^ m + Oih 2 ) . 

Special attention will be given to the smoothness of the Oih 2 ) coefficients. 

To show (6.8b) we start by using the trapezoidal rule to approximate the 
integral in (2.1b); we get 

(6.9a) fj+j2( u ) = yl/(“( x ;+i/2> f n)) + /(“(*/+ in* + t))] + Oih 2 ) . 
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The smoothness of the O^h 2 ) term follows from that of /(a) and u(x,f). Next 
we observe that a(x) in (6 4b) approximates a(u(x,f„ )) to 0(h?), and there¬ 
fore we can use the approximate characteristic line (6.5c) to trace 
u(x ; + 1 / 2 , t n + t) to u(x Q (x j+l/2 , t), t n ) with 0{h 2 ) accuracy; consequently 

(6.9b) f(u(xj+v2, t n + t)) = /(«(x 0 (x J+Wl t),' t n )) + 0{h 2 ) . 

Finally we obtain (6.8b) by approximating u(x,t n ) in (6.9a) and (6.9b) to 
0(h2) by L(x; u n ) (see (4.4)). The smoothness of the Oiji 2 ) term in this 
approximation is due to (4.4c): 

SJ = A • «,(*, (.) + Oik 3 ) . 

(We recall that the degeneracy to first order accuracy at local extrema points of 
some "second-order accurate" TVD schemes is due to lack of smoothness there of 
the O^h 2 ) term in (2.7a)). 

We turn now to verify (6.8c). First let us consider the case aj+y 2 ^ 0: 

£(*>+ 1 / 2 ; “") = “7 + \ SJ . 


Lifrfoj+yb r); O = “? + 


\a 


y+ 1/2 


1 + X(fljM.i/2 - *jf-l/2) 


and expand the LHS of (6.8c) around uj. We get 

(6.10a) /( uf) + | <x(57)(l - X «, +w ) S, + ja'(5fl[(l - X a,.^ 2 


+ + 0(A 3 ) - /j+yj + Y (Zk 2 a 2 - 1) ■ a' • (a,)V l/z + 0(A 3 ) . 


Similarly in the case s 0: 
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*-(*/+!/ 2i •*")*“/+1 - i -SJ+1. £(*c(*)+1/2> t ); 5") 


— iT 71 

~ «V+i 


* fl y+i/2 


2 1 + M fl /+3/2 - ®jf*l/2) 




we expand the LHS of (6.8c) around 5jf +1 to get 
(6.10b)/(uy +1 ) - y a(tty +1 )(l + kaj+y^Sj+i 

+ j ®'(“7+i)[(l + x «/+m) 2 + ( X ^*- J i) 2 ] • {Sj+ 1) 2 + 0(A 3 ) 


- //+w + Y (2AV - 1) • a' • WW + 0(* 3 ) . 

We see from (6.10a) and (6.10b) that independently of the sign of a ;+1/2 » the 
0(A 2 ) term in (6.8c) is the same, namely 

Y ( 2x2flZ “ x ) * ^ * (“x) 2 |/+i/2 • 

This completes the proof that the scheme (6.7) is second-order accurate in the 
sense of (2.9) wherever u(x,t) is smooth, including local e x tr em a and sonic 
(f = 0) points. 

Remarks: (1) The numerical flux (6.7b) can be rewritten as: 

(6.11) f J+u2 = |[/(v7) +/(vf +1 ) - F J+W I (v7 +1 - vf) 


+ max(0, ay + y 2 ) ' (1 “ 2 ) ' $j 


— min(0, 5/+i^) * (1 + * 
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(2) Our proof that the scheme (6.7) is non-osdllatory is based on the 
representation of (6.7) as the cell-average (6.2) of the non-osdllatory approxi¬ 
mate solution v„(x,t) in (6.6). To ensure that v„(x,i) remains univalued for 
OifSr we have to restrict the time-step t so that for all j 

(6.12a) */+i/ 2 (t) > xj-wi t) . 

This implies the condition 

(6.12b) t max j (a j ^y 2 “ S/+ 1 / 2 ) ^ h • 

On the other hand, to derive the particular form of the numerical flux (6.7b) we 
have made the assumption 

(6.13a) \xj +i/ 2 (t) - xj+vi I < h for all ; , 

which implies the condition. 

(6.13b) t max^lay*< h . 

The two time-step restrictions (6.12b) and (6.13b) are characteristic to 
Godunov-type schemes. The practice of reconciling the two conditions by 

t maxy|a]i +1 /2| ^ h /2 

is completely unnecessary: The scheme (6.7), as the original Godunov scheme, 
remains non-osdllatory (or TVD in the case analysed in [10]) for the full CFL- 
restriction (6.13b). The reasoning for this observation is as follows: The approxi¬ 
mate solution (6.6) under the CFL restriction (6.13b) may fail to be univalued in 
they-th cell only if a)- 1/2 > 0 and a/+i /2 < 0. In this case the numerical fluxes 
fj ± i /2 as defined by cell-averaging in the neighboring cells j ± 1, remain (6.7b). 
Thus the only thing that needs to be changed in this case is the definition of 


v„(x,t) in the j-th cell. 

(3) We observe that once v„(x.t) is defined globally in (6.6) there is no 
intrinsic need to average it on the original mesh. We may average it on different 
intervals and still conclude that the resulting approximation is non-osdllatory and 
conservative. Furthermore, the construction of the interpolant Q, the approxi¬ 
mation L and the approximate characteristic field a(x) needed to define 

v„(x, t), does not depend on the uniformity of the mesh. Therefore the scheme 
(6.7) generalizes immediately to non-uniform moving meshes. Of particular com¬ 
putational interest are the self-adjusting moving grids of the type described in 
[12], which make it possible to obtain perfectly resolved shocks and contact 
discontinuities. 

(4) We note that since the approximate solution v„(x,r) in (6.6) is conser¬ 
vative, it is possible to consider an associated random-choice method obtained by 
replacing the cell-averaging in (6.2) by a sampling with respect to a variable that 
is uniformly distributed in the cell, i.e. 

V/ +1 - V„(x y + 8 /A,t) 

where Qj is uniformly distributed in [-1/2,1/2]. Following the reasoning of [7] 
it is clear that the resulting random-choice method is non-osdllatory and that its 
limi ts are weak solutions of ("Li). Although the random-choice approach has 
many attractive computational features, it has been our experience that in many 
application it is possible to accomplish the same computational goals with a self- 
adjusting moving grid. In this case the use of the latter is preferable as it offers 
gain in resolution without a loss in pointwise accuracy that is associated wL\ sam¬ 
pling. 
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7. Numerical Illustration. In this section we compare the new uniformly 
second-order non-ounllatory scheme of this paper (to be refered to as UN02) to 
the typical second-order TVD scheme (to be refered to as TVD2). Both schemes 
can be written in the form (6.7), i.e. 

(7.1a) vj+ l • vj - ^(fj+vi - /j-\rd » 


(7. lb) 1/2 


/(vj 1 ) + 1/2 5)+ 1 / 2(1 ~ + X(a^4.i/2 “ «/- 1 / 2 )] 

if ^f+i/2 ^ 0 

/O^+i) ~ 1/2 ^+ 1 / 2(1 + Vi ~ ajf+i/2)l 

if ^jf-M/2 s 0 


(7.ic) sj - 

here 5) +1/2 is (6.4a) and m(x,y) is the min mod function (3.3). Sf are dif¬ 
ferent for TVD2 and UN02: 

(7.2) TVD2: Sf ^d J±u2 v n 


(7.3) UN02: Sf = ^ lU2 v” T y D J±yl v" , 

where and are defined in (3.2). 

UN02 and TVD2 are both second-order accurate Godunov-type schemes 
that differ only in the reconstruction step (4.2a): 

(7.4a) L(xyi) = u< f Sj(x - xj)/h for \x - xj\ < h/2 , 

when 1 ! the slopes of the lines are calculated by (7.3) and (7.2), respectively. 
Therefore we start our comparison on the approximation level. 

In Table 1 and Fig. 1 we present approximations to 
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u(x) * simrx, -lsxsl. Wc divide [-1,1] into V equal intervals and 
define 

(7.4b) ZJ--1+J-1. 0 sjxN. 

The symbols in Fig. 1 denotes values of uj - sinux^ for N * 10 in (7.4b). In 
Fig. la we show the piecewise-parabolic interpolant Q(x\u) (see section 3). In 
Fig. lb we show the piecewise-linear approximation L^x^) which is (7.4a) 
with (7.1c) and (7.3). In Fig. lc we show the piecewise-linear approximation 
L TVD2 (x;u) which is (7.4a) with (7.1c) and (7.2). We make the following obser¬ 
vations regarding Fig. 1: (i) Q is a better approximation than I UN02 ; I UND2 is a 
better approximation than L TVD2 . (ii) rV(L UN02 ) > 7V(u) > 7V(I TVD2 ). hi 
table 1 we quantify the first observation; we list the L«-error and the 12 -error of 
these approximations to simrx for a refinement sequence of N = 10,20,40,80 
in (7.4b). Gearly Q is an 0(h 3 ) approximation, while L UN02 and L TVD2 are 
0(/ 1 2 ). The error in L UN02 is about a 1/3 of the error in L TVD2 . 

In Table 2 and Fig. 2 we present solutions of UN02 and TVD2 for the con¬ 
stant coefficient case 

(7.6) u, + m, = 0, u(x,0) = sin irx, -1 as x ^ 1 

with periodic boundary conditions, at t = 2 with r/h - 0.8. Figs 3a and 3b 
show UN02 and TVD2, respectively, with N * 20 in (7.4b). In table 2 we list 
the Loo-error and Lj-error for a refinement sequence with N - 20,40,80,160. 
Gearly UN02 is second-order accurate in both L m and L lt while TVD2 is 
second-order accurate in L\ but only first order accurate in L«. Fig. 3b demon¬ 
strates that the degeneracy to first order accuracy at local extrema in the TVD 
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scheme adversely affects the accuracy everywhere (Because the scheme is TVD it 
cannot switch ab. aptly to second-orr*T accuracy as this would introduce osdlla- 
tious; consequently it takes quite a while to recover the second order accuracy). . 

Next we approximate the discontinuous function 

* 

-xsin(3irx 2 /2) , -1 < x < -1/3 
(7.7) u(x) - |sin(2'irx)| , |x| < 1/3 

2x - 1 - l/6sin(3irx) , 1/3 < x < 1 

which we extend to have period 2 outside [-1,11- 

In Fig 3 we present approximations to <J>(x), using N = 20. Fig 3a shows 
C(x;u), Fig 3b shows L UND2 (x;5), and Fig 3c shows L TVD2 (x;u). We again 
observe that Q is better approximation then L UN02 , while L lJS02 is a better 
approximation then L TVD2 . 

In Fig 4 we present solutions of UN02 and TVD2 for the constant coeffi¬ 
cient problem (7.6), initial data given by (7.2), and periodic boundary conditions. 
We take t = 2 and j/h =* 0.8. Figs 4a and 4b show UN02 and TVD2 respec¬ 
tively with N = 40. Fig 4b shows the damping effect that the TVD scheme 
imposes due to its degeneracy to first order accuracy at local extrema. 

In Fig 5 we solve the same problems, except we impose boundary conditions. 
At x = -1 we impose the given function (7.7) evaluated at -1 - r. No boun¬ 
dary conditions are imposed at x = 1. We implement this numerically using 
UN02 and TVD2 except at the boundary points. There we are in general, unable 
to construct non-osdllatory piecewise parabolic interpolants Q(x,u), so we con¬ 
struct the only possible parabolic intr.polant thru x it Xj+i and the point to either 
the left or right which lies in the region. The analogous procedure is carried out 
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at the reconstruction stage. Figs 5a and 5b again show the results at t - 2 with 
rlh = 0.8. 


The possible introduction of oscillations through the boundary conditions 
does not seem to have degraded the p e rformance of either scheme (in fact the 
opposite is observed). Again the TVD2 scheme shows a damping effect. 

In Table 3 and Fig. 6 we present results for Burgers’ equation 


(7.8) u, + uiij - 0, u(x,0) - a + sin n^x + 0), -lsxs 1 


with periodic boundary conditions and r/h (1 + |a|) - 0.5. The solution to (7.7) 
is smooth for t < 1/u; at t * 1/ir it develops shocks. In Table 3a we list the 
L.-error and Lj-error of UN02 and TVD2 at t - 0.15 for a = 0 = 0 in 
(7.7). This table shows the same asymptotic behaviour as Table 2, except that 
because of the large gradients it shows for a smaller h. 

In Figs 6a and 6b we show results of UN02 and TVD2 for (7.8) with a = 2 
and 0 » 1 at / « 0.35 with N =* 20. In this case the solution to (7.8) 
develops a shock moving with speed 2 beginning at time r » 1/ir = 0.318. 

In Table 3b and Figs 6c and 6d we repeat the previous calculations for the 
schemes (2.13): 

(7.9a) vf +1 = v7-X( fj+M-fj-vd 


(7.9b) fj+y 2 * /(v„(x ;4 . V2, t/2)) =■ f(L (xg(xy4 .1/7, t/ 2 ),v h )) 

where t/2) is (6.5c), i.e. 


VWj+vi = 



J_ 1 - s „) 

2 1 + XOfy+i/z - dj. in y2 J J 

1_ 1 + H&j+yi + d)+\nV2 

2 1 + \ 


if 1/2 ^ 0 

if 4/* 1/2 ^ ® 
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As we have remarked in section 2. v}* 1 in (7.9a) is not a cell-average of 
v„(jr,r), but only an approximation to it. Therefor* it is not necessary to take 
dj+v2 (7.8c; to be (6.4a). We choose so that (7.9c) is continuous at 

&)+vi * 0 : 


(7.9d)J^ w 


- 7-%)-/(>•;+ 


(v?*i - - (vj 1 + y57) 


We denote the schemes (7.9) with SJ defined by (7.1c) and either (7.2) or 
(7.3) by FVD2 and FN02. respectively. We note that (7.9) is identical to (7.1) in 
the constant coefficient crse, and consequently FVD2 and FN02 are nonosdlla- 
tory in the constant coefficient case. Figs 6c and 6d show that FN02 and FVD2 
are also non oscillatory in the case (7.8). Furthermore, Table 3b shows that 
FN02 is much more accurate than UN02 (FVD2 is about the same as TVD2). 

In all previous examples we have pr es ented pointwise calculations; namely, 
we have initialized the numerical solution by taking vP to be the value of the ini¬ 
tial data at x ; , and we have considered vj to be an approximation to u(xj,t„). 
(Surely this is an acceptable practice for second order accurate schemes.) In Table 
3c we repeat the calculation for UN02 in Table 3a, but now in a sense of cell- 
averages and denote it by AN02. Now we initialize UN02 for (7.8) with 
a = 3 = 0 by cell-averages of the initial data, i.e. 

(?.10*)vf «« » (cos(ir xj+yfi - cos^Xy-^] = x y ) , 

and regard vj to re p resent cell-averages of u(x,f„). To obtain a pointwise 
approximation to u(x,r„) we first compute point values ^ of its indefinite 

integral u n {xj^ =■ f ,t H ) dy by 

Xq 
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(7.10b) 


V7*w-A f vf- 

l-io 

Next we obtain a global piecewise-linear approximation v(x,/„) to u(x f t„) by 
(7.10c) »(*,/„) = ■£ 

where Q is the piecewise-parabolic interpolant of section 3. Finally we get 
(7.10d) v(xj,t.) - Q(xjin = j (Vf +U2 - '/;- w ) = yf . 

Tlius the only difference between AN02 in Table 3c, and UN02 in Table 3a 
is the initialization (7.10a), which itself differs only slightly from the mesh values 
of the initial data (since sin(irA/2)/(irA/2) = 1 - l/6( , rrA/2) 2 + 0(h 4 )). 

We remark that cell-averages do play a significant role when the initial data 
is discontinuous (since they provide information about tne location of the discon¬ 
tinuity) and in higher-order Godunov-type schemes; this will be described else¬ 
where. 

We turn now to present calculations with a formal extension of UN02 and 
TVD2 for systems of conservation laws. We consider a Riemann problem for the 
Euler equations of gas dynamics 

( tt£ x < 0 

Ur x >0 » 

(7.11b) u = (p,m,£) r ,/(u) = (m,m 2 /p + P, m(E + ?)/p) r , 

(7.11c) P = (y - 1 )(E ~ y m 2 /p) . 

Here p, m, E and P are the density, momentum, total energy and pressure, 
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respectively; we take y = 1.4. 

In the following we apply the extension technique of [3] to UN02 and TVD2. 
The idea is to extend UN02 and TVD2 to systems in such a way that will be 
identical to (7.1) in the scalar case, and will decouple into (5.3) for each of the 
characteristic variables in the constant coefficient system case. To accomplish 
that we use Roe’s averaging for (7.11) (see [13]) 

(7.12a) Vj+w = V(vJ, v/ +1 ) 

for which 

(7.12b) /(yj+i) “ f(yj) = a ( v 7+i^)( v 7+i “ v 7). A(u) = 3//3u , 

and define local characteristic variables with respect to the right-eigenvector sys¬ 
tem of A(vy +1/2 )- We extend (6.11) to systems as follows: 

(7.13a) v/ +1 = v;->4 +1/2 - fj. i/j) 

(7.13b) j=| j/(»f) + ^ cf +yl Rf +in 

(7.13c) Cj+i /2 = \aj+\/2\df+V2 w ~ max (0,fly : + i/2)(l — kaj- in)Sj 

+ min(0,ot + i/2)(l + ^^+3/2)*^/+i • 

Here aj+ \n. is the *-th eigenvalue of <4 (vj + 1/2 ) corresponding to Rf+uz, and 
dj+ \a w denotes the component of dj+y$y = v™ +1 - v" in the k-th characteris¬ 
tic field, i.e. 

(7.13d) dj+iftV = 2 (dj+y2 w )Rj+\n • 
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Likewise Sj denotes the component of the vector of slopes in the ifc-th charac¬ 
teristic field, and is defined as follows: 

(7.13e) Sj = m(5* j,5*j)/[ 1 + k(a ^ +1 /2 ~ ^- 1 / 2 )] ] 

m(x,y) is the min mod function (3.3). S±j are different for VD2 and N02: 

(7.14) TVDl:S k ±J = d} ±y2 w 

(7.15) UN02: S k ±J = d} +vl w =F ± Df ±ul w, 

Di±U2 w = *(<$+3/2 W “ 1/2*. - <#-1/2*0 • 

In Figs 7.8 and 7.9 we show numerical solutions of UN02 and TVD2, respec¬ 
tively, for the Riemann problem (7.1b) with 

U L = (l,0,2.5) r , U R = (0.125,0,0.25) . 

These figures demonstrate that the formal' extension to systems is nonoscillatory 
in this case. Since the solution to the Ri emann problem is just constant states 
seperated by waves we do not get to see here the extra resolution power of UN02, 
except that its numerical solution is somewhat "crisper" than that of TVD2. In 
this calculation we have not employed any artificial compression in the linearly 
degenerate field and therefore the contact discontinuity smears like n 173 , as 
expected. The interested reader is referred to [4], [5] and [10] for a detailed 
description of such compr es sion techniques, as well as for details of entropy 
enforcement mechanisms. 
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Table 1: Approximations to u(x) = sin(‘irx), -lsisl, with periodic boundary conditions. 



Loo-ERROR 

L r ERROR 

N 

Q 

^UN02 

^TVD2 

Q 

^UN02 

L TVTJ2 

10 

| 

B * 

1.420 x 10” 1 



|^TOTW| 

20 

1.971 x 10 -3 

1.231 x 10 -2 

3.558 x 10" 2 

1.802 x 10~ 3 

5.576 x 10" 3 

1.525 x 10" 2 

40 

2.476 x 10" 4 

3.083 x 10“ 3 

9.163 x 10~ 3 

2.148 x 10“ 4 

1.355 x 10" 3 

3.902 x 10 -3 

1 

3.104 x 10“ 3 


nimri 

2.617 x 10“ 3 

3.351 x 10" 4 

9.787 x 10“ 4 


Table 2: Numerical solutions of u, + u, = 0, u(x,0) = sin wx, -1 s x s 1 at t = 2 
with periodic boundary conditions and r/h = 0.8. 



L.-ERROR 

Li-ERROR 

N 

UN02 

TVD2 

' UN02 

TVD2 

20 

7.097 X 10" 3 

8.119 x 10 -2 

8.944 x 10" 3 

6.778 x 10“ 2 

40 

1.607 x 10" 3 

3.477 x 10“ 2 

2.044 x 10“ 3 

2.033 x 10" 2 

80 

3.870 x 10" 4 

1.453 x 10“ 2 

4.926 x 10" 4 

5.626 x 10“ 3 

160 

9.201 x 10” 5 

5.975 x 10" 3 

1.172 x 10“ 4 

1.528 x 10“ 3 
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Table 3a. Numerical solutions of u, + uu x = 0, u(x,0) = simrx, at t = 0.15 and r/h ~ 0.5 - 
UN02 and TVD2 



L<»-ERROR 

L r ERROR 

N 

UN02 

TVD2 

UN02 

TVD2 

20 

", 


| 


40 

5.712 x 10" 3 

1.054 x 10" 2 


5.051 x 10’ 3 

H 

1.552 x 10“ 3 

4.422 x 10" 3 


1.340 x 10" 3 

160 

3.985 x 10” 4 

1.837 x 10“ 3 

1.965 x 10“ 4 

3.621 x 10” 4 


Table 3b. Same as Table 3a for FN02 and FVD2. 



L«-ERROR 

Lj-ERROR 

N 

FN02 

FVD2 

FN02 

FVD2 

20 

| 


| 


40 

1.959 x 10" 3 

1.054 x 10" 2 

8.869 x 10"“ 


m 

5.106 x 10“ 4 

4.424 x 10" 3 

2.163 x 10~ 4 

1.072 x 10" 3 

160 

1.251 x 10" 4 

1 


2.946 x 10" 4 
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PIECEWISE PARABOLIC 



Figure la 

~i-rure 1 . Approximations of u = sinnx. -1 ^ x * 1, 
a) Piecewise-Parabolic interpolant Q(\,u). 


with N - 10. 



Plecewlse-llnear approximation L TVD2 



PIECEWISE LINERR TVD PIECEWISE LINERR 




















ngureJia. Numerical eolation of u,. + U)[ - 0, u(x,0) Figure lib . Same as Figure la using TVD2. 

defined by (7.7) at t -« 2 with periodic 
boundary conditions N = *40 t/h - 0.8, 
uainR UN02. 











Figure 6a. Solution to (7*8) wiih o 
CFL if =.y, at t = 0.35 using UND2. 






Figure 6C, Same as Fig. 6a using FN02. Figure 6P . Same as Fig. 6a using FVD2. 





Figure 7b . 

Numerical solution of density in a Riemann problem for Euler's equations, 
(b) TVD2 















fijirerical solution of velocity in a Rlemann problem for Filler's equations. (a) UIJ02. (b) T D2. 
















